!==========================================
! 28.11.12, LNGS
! rewritten for FORTRAN 90
! 
module mlt
   implicit none

!------------------------------
     contains
!------------------------------


SUBROUTINE MULTY(EXT,RES,NN,M,N,EPS0)
   implicit none
   integer,intent(in):: M, N
   integer,intent(in),dimension(N):: NN
   real(8),dimension(M):: RES
   real(8),intent(in):: EPS0
   real(8):: EPS, UK, AA, BB, DW
 
   INTEGER P, Q
   integer,dimension(4):: H, K, S  
   integer:: L, I, I1, J, JJ, K1, K2, E, E1

   real(8),DIMENSION( 4):: A, D, X, U, W
   real(8),DIMENSION(10):: Y, V
   real(8),DIMENSION(4,10):: F, R
   
   DATA U/0.8611363d0, 0.3399810d0, 0d0, 0d0/   &
      , W/0.3478548d0, 0.6521452d0, 0d0, 0d0/   &
      , P/4/
!- - - - - - - - - - - -
   EPS=EPS0
   DO JJ=1,N  
      S(JJ)=MAX0(1,NN(JJ))
   enddo  ! JJ

   DO I=1,M
      V(I)=0d0
   enddo
   
   L=0
   DO JJ=1,N
14    IF(L.EQ.1) S(JJ)=2*S(JJ)
      L = 1
      J = 0
      GOTO 3
!- - - - - - - - - - - -
1     K(J) = 1
2     K1 = K(J)
      K2 = P+1-K1
      UK = -U(K1)
      IF(K1.GT.K2) UK = U(K2)
      X(J) = A(J) + (2d0*H(J) - 1d0 + UK)*D(J)
      IF(J.EQ.N) GOTO 5
3     Q = J+1
      DO J=Q,N
         CALL EXT(J,X,AA,BB,Y)
         A(J) = AA
         D(J) = (BB-AA)/(2*S(N+1-J))
         X(J) = AA + (1d0 - U(1))*D(J)
         H(J) = 1
         K(J) = 1
         DO I=1,M
            R(J,I) = 0d0
            F(J,I) = Y(I)
         enddo ! I
      enddo ! J  
      J = N
5     CALL EXT(N+1,X,AA,BB,Y)
      K1 = MIN0(K(J),P +1 - K(J))
      DO I=1,M
         R(N,I) = R(N,I) + Y(I)*W(K1)
      enddo ! I
7     K(J) = K(J) + 1
      IF (K(J).LE.P) GOTO 2
      H(J) = H(J) + 1
      IF(H(J).LE.S(N+1-J)) GOTO 1
      J=J-1
      IF(J.EQ.0) GOTO 9
      K1=MIN0(K(J),P+1-K(J))
      DW = D(J+1)*W(K1)
      DO I=1,M
         R(J,I) = R(J,I) + R(J+1,I)*F(J+1,I)*DW
      enddo ! I
      GOTO 7
9     DO I=1,M
         Y  (I) = S(I)
         RES(I) = R(1,I)*F(1,I)*D(1)
      enddo ! I
!     - - - - - - - - - - - -
      CALL EXT(N+2,Y,EPS,EPS,RES)
      IF(EPS.LE.0d0) RETURN
      DO I=1,M
         E  = DABS(RES(I)-V(I))
         E1 = DABS(EPS0*RES(I))
         IF(E.LE.E1) EXIT
         DO I1=1,M
            V(I1) = RES(I1)
         enddo ! I1
         GOTO 14
      enddo ! I
      S(JJ) = S(JJ)/2
   enddo ! JJ
   RETURN
END SUBROUTINE MULTY


end module mlt
!++++++++++++++++++++++++++++++++++++++++++++++++

!----------------------------------------
Program tstmlt
   use gazlib, only: multy, dgauss
   implicit none

   integer,dimension(1):: NN
   external FFF
   real(8),dimension(2):: RES
   real(8):: EPS0 = 1d-8, rs

   NN(1) = 0

   call MULTY(FFF,RES,NN,2,1,EPS0)
   res
   write(*,'(2E12.5)') RES

   rs = dgauss(Fsn,0d0,1d1,1d-8)
   write(*,'(E12.5)') rs

end Program tstmlt

!----------------------------------------
Subroutine FFF(J,XX,XL,XU,R)
   implicit none
   integer:: j
   real(8),dimension(1):: XX, XL, XU
   real(8),dimension(2):: R
   real(8):: X
   goto (1,2,3), j
1  XL(1) = 0d0
   XU(1) = 1d1
   R(1)  = 1d0
   R(2)  = 1d0
   return
2  X = XX(1)
   R(1) = x*sin(x)
   R(2) = x**3
   return
3  return
end Subroutine FFF

!----------------------------------------

real(8) Function Fsn(x)
   implicit none
   real(8),intent(in):: x

   Fsn = x*sin(x)

end Function Fsn

